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Abstract 

In this paper the revised Kajantie-Byckling approach and improved phase space sampling techniques for the massive 
multi-particle final states are presented. The application of the developed procedures to the processes representative for 
LHC physics indicates the possibility of a substantial simplification of multi-particle phase space sampling while retaining a 
respectable weight variance reduction and unweighing efficiencies in the event generation process. 
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1 Introduction 

With the advent of LHC era a need for precise predictions (and subsequently detailed simulation) of many QCD and elec- 
troweak processes has arisen. One of the necessary components of an accomplished simulation tool is definitely efficient 
phase space modeling of multi-particle final states. While the light (massless) particles in the final states and the relevant 
topologies can quite effectively be simulated using general tools (e.g. RAMBO,SARGE,HAAG 0|4]|2l) further develop- 
ment might be needed for phase space sampling where the massive final state particles are present. While at LEP's 200 GeV 
a massless approximation for final state particles was adequate, in most of the studied cases at LHC in contrast the massless 
approximation becomes less suitable because of crossing the top quark production barrier in conjunction with the possibility 
of multi-quark final states and the shifting centre-of-mass energy in proton-proton collisions. Furthermore, in the LEP era the 
number of (hard process) particles in the final states of relevance only rarely rose above four and subsequently the number of 
possible final state topologies was relatively modest. Typical processes of interest (signal) and their backgrounds at LHC will 
involve several heavy quarks and/or weak bosons in the intermediate states; thus the number of particles in a typical process 
under study is generally at least four. Furthermore, the number of Feynman diagrams contributing to a typical process involv- 
ing heavy quarks in the final state is mostly ranging from a few to hundred(s) (and can steeply rise to many thousands when 
massless quarks are added). Subsequently, this results in varied particle topologies which prove to be a challenge when trying 

*Partly supported by Marie Curie Host Fellowship for the Transfer of Knowledge Contract No. MTKD-CT-2004-510126 and by the EC FP5 Centre of Excellence "COPIRA" 
under the contract No. IST-2001-37259. 
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to adequately describe them using the available statistical approaches and numerical methods. Important steps have already 
been made in matrix element calculations (e.g. MadGraph|6 1) and have so far surpassed the corresponding development of 
phase-space modeling and sampling techniques. 

The general objective in simulation of physics processes for the LHC environment is thus to improve the integration of 
the differential cross-section using Monte-Carlo sampling methods 1 . The sampling method used should aim to minimise the 
variance of the integral as well as maximise the sampling efficiency given a certain number of iterations and the construction of 
the sampling method itself should aim to be sufficiently general and/or modular to be applicable to a wide range of processes. 
Writing down a (process) cross-section integral for LHC type (hadron-hadron) collisions: 



a = l52 Mxi,Q 2 )f b (x2,( 



\Mn\ 



(2tt) 3 "- 4 (2s) 



d,Xl dX2 d$ri, 



(1) 



where f a ,t>(x, Q 2 ) represent the gluon or (anti)quark parton density functions, |-M n | 2 the squared n-particle matrix element 
divided by the flux factor [(2-7r) 3n_4 2s] and d$ n denotes the n-particle phase space differential. The quantity s = xi X2 s is 
the effective centre-of-mass energy, and the sum X] a b mns m case °f q uar k- ant iquark incident partons over all possible quark- 
antiquark combinations (a, b = u, d, s, c, u, d, s, c). In case of gg initial state the sum has only one term with a = b = g. 
It is often convenient to re-write the differential cross-section in the form: 



Xl, 



") X' 2 fb(x2, Q ) 



\Mn\ 



(27r) 3 "- 4 (2s 2 ) 



dy ds 



(2) 



with the new (rapidity) variable given by y = 0.5 log(xi/x2). The n-body phase-space differential d<E> n and its integral 
$n depend only on s and particle masses m, due to Lorentz invariance: 

$n(s,mi,m 2 , . . . ,m„) = J d& n (s,m 1 ,m 2 , . . . ,m n ) — J S 4, ^(p a + p b ) - ^P'^ ]~[ d i p i S(p 2 - m 2 )9(p°), (3) 

with a and b denoting the incident particles and i running over all outgoing particles i = 1, . . . ,n. What one would like to do 
is to split the n-body phase parameterised by 3n-4 essential (i.e. non-trivial) independent variables into manageable subsets 
(modules) to be handled by techniques which reduce the variance of the result and/or the sampling efficiency (e.g. importance 
sampling|7 1 or adaptive integration like VEGAS 1 8| or FOAM| 9 1). Stating this in formal terms, the above Equation|2|should 
be transformed into an expression like: 



/ 



n / ds * 



\ 
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where one integrates over Mandelstam type (Lorentz invariant) momenta transfers Si, tj and space angles Qk = (cos t?fc, (j>k) 
within the kinematically allowed limits (3n-4 variables in total) with the term \J n \ denoting the Jacobian of the transformation. 
If one would then decide to introduce importance sampling functions in order to reduce the peaking behavior of the integrand 
1 7 j, the integrals would take the form: 



dsi = 



9i(si) 



dsi, 



(5) 



where the importance sampling function gi is probability density function normalised in the integration region [s ; , g; 

h 

gi{si)dsi = 1, 



(6) 



'For a nice discussion on the topic see e.g. 1T1121 



B. P. Kersevan,E. Richter-Was: 



Improved Phase Space Treatment . 



3 



which exhibits a similar peaking behavior as the integrand. Formally, one then inserts the identity: 

i f H \ 



gi(si)dsi 



\ 



dn 



(7) 



/ 



into the integral and then derives the unitary sampling prescription: 

- / gi{si)dsi 



i 



dr, 



(i 
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/ gi(si)dsi — - — r-dsi = [ dri [ 5 (si — G (r»)) — 7 - 
J 9i{Si) J J v ' gi(s-. 

/ 



-dsi = 



dri 



(8) 



/ 



which formally means that the Si values are sampled from the interval according to the gi (si) distribution by using (pseudo-)random 
variable r; together with the gi(si) cumulant G(si) = J s l gi(si)dsi with its inverse G -1 . The unitarity of the algorithm states 

that each trial ( n value) produces a result (i.e. a corresponding sj value distributed according to gj(sj)). 
Performing such substitutions on all integration parameters would give as the cross-section expression; 

/(n,r 2 . . .) 



a = IT I dr 



' g(ri,r 2 ■ ■ .) 



(9) 



where the integrand would (hopefully) have as low variation as possible at least for a subset of contributing Feynman dia- 
grams 2 . To improve the sampling method further, the n (pseudo-)random variables can be sampled from adaptive algorithms 
of the VEGAS type ISlfTol. 




Fig. 1: A representative Feynman diagram describing a 2 — > 6 process gg — *• tt — > bbW + W 
bblviq\q2 and its decomposition into a set of 2 — + 2 t-channel and s-channel sub-processes. 



A representative Feynman diagram describing a 2 — » 6 process is shown in FigureQ] As one can see, the process can be 
split in several consecutive branchings, this approximation is often used in matrix element (probability amplitude) calculations. 
It seems rather obvious that any Feynman diagram can be split in a series of horizontal and vertical branchings that one can 
denote as s-type and t-type(u-type) using the analogy with the Mandelstam variables. What one would like to do is thus to 
modularise the phase space in the form of sequential s- and t-type splits. 

The s-splitting of phase space is relatively easy to do and has as such been used in many instances of Monte-Carlo 
generation (e.g. FermiSV 1 12 1, Excalibur 1 13 , Tauola 1 14 1 etc.); the t-type branchings (often tagged as multi-(peri)pheral 
topologies) have in contrast generally been calculated only for specific cases (e.g. for 3 or 4 particles in the final state 1 12l ll3l ). 
As it turns out, the problem of several massive particles in the final state has already appeared more than 30 years ago when 
several hadrons (e.g. pions) have been produced in (comparatively low energy) nuclear interactions. At that time Kajantie 
and Byckling 1 15 1 have derived the formulae for simulating any sequence of s- and t- type branchings which, with some 
modifications, can also be applied to the EW and QCD processes involving heavy quarks and/or massive bosons at LHC. 

In the following Sections |2| and [5] the revised version of of Kajantie-Byckling algorithm (KB) will be discussed and a 
formulation of the algorithm for the multi-phase space integration will be presented. In Section |4] the numerical results for 
some representative applications as implemented in the AcerMC 2.0 Monte-Carlo generator will be presented. 



The 'modularisation' can be performed for several topologies at the same time and multi-channel techniques can be applied. 
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2 Modified Kajantie-Byckling Formalism 

2.1 The s-type Branching Algorithms 




Fig. 2: A diagram of a generic 2^2 s-channel process. 

The s-splits are the simplest method in the KB formalism. For the sake of completeness one should start with the definition 
of the two-body phase space integral (c.f. Fig[2j: 

$ 2 (s,rai,ra 2 ) = J d 4 pid 4 p 2 5(p 2 1 - ml)5(pl- ml)S 4 (p- p! - p 2 )0{pi)0{p2), (10) 

with the incoming momentum sum p = (p a + pb),p 2 = s and outgoing momenta pi,2,p 2 ,2 = m i ,2- The phase space 
integral is Lorentz invariant (as one can observe in the above Equation where it is written in a manifestly Lorentz invariant 
form). Subsequently, due to Lorentz invariance, the integral is necessarily a function of the Lorentz scalars s, mi and ni2 only. 
The step function product 0{pi)G>(p (1 ) is the explicit requirement of the positiveness of the energy terms in pi,2 while the 
delta functions represent the on-shell conditions on pi,2 and the total momentum conservation. 

The integral can be transformed into a more compact form by integrating out the spurious variables; one thus first integrates 
over d 4 p2 and chooses the centre-of-mass system (CMS) as the integration system of reference with p = (*/s, 0, 0, 0) and 
then evaluates the integrals over p? and Ej : 



<£>2(s,mi,m 2 ) = J d A p 1 8(p\ - m\)8((p - pi) 2 - m|)9(p?) 



(11) 



.1-1 ft."--"! I J^l „ /— 

pt(s,mi,m 2 ) 



with the stars explicitly denoting the values in the centre-of mass system. The first integration simply sets p? = (Pi) 2 + m 2 = 
E* and the second integral leads to the well known relations for the energy: 

, _ s + ml - ml , _ r , _ s + m\ - mj 

and momenta sizes: 



Pl = \Pl\ = , P2=Pl (13) 

of two particle production. The A(s, m 2 , m|) denotes the Lorentz invariant function: 

A(s,m 2 , ml) = (s — (mi + m 2 ) 2 )(s — (mi — m 2 ) 2 ) (14) 

and thus explicitly contains the phase space cutoff, i.e. the requirement that the available CMS energy y/s should be bigger 
than the mass sum y^s > (mi + 1112). Note that the integration was so far done only over the spurious parameters, leaving the 
polar and azimuthal angle of the pi particle as the two independent parameters d£T = dcos0*d(p* . The integral becomes 
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trivial to sample in case the outgoing particles can be approximated as massless (the 'boost' factor lambda transforms to 
unity). As already claimed, the latter approximation is however often unjustified when studying processes representative for 
the LHC environment. 

Kajantie and Byckling 1 15 | introduced the recursion and splitting relations for the n-particle phase space $ n (s) given by 
Eq.|3] The recursion relation can be derived by defining the momentum sum: 



hi — ^ Pj — (hi , ki)] Mi — ki 
3=1 



Subsequently one can interpret p = k n and s — M„ from Eq.|3| One continues by introducing the identities: 

i = J dM^sikl-! - M*_i)e(j£_ x ) 



and 



1= / d 4 k n -i8 4 (p — k 



n-l — Pn) 



(15) 

(16) 
(17) 



(18) 



into the integral of Equation^ separating out the arguments containing k n _i and p n terms one obtains: 

$„(A^,mi,m 2 ,...,m„) = J dMn-i X 

x I J d i k n - 1 d i p n 5(k' 2 n _ 1 - Afn_i)5(p* - ml)8 4 (p- k n -i -p«)6(A;°_i)6(p°)| x 

x $„_i(M^_ 1 ,mi,m 2 ,...,m„_ 1 ), 

where the remaining pi terms form the (n-l)-particle phase space integral $ n -i (M^_i, mi , m2, . . . , m n -i) and the terms in 
curly brackets give a two particle phase space term (c.f. Eq. ll2> : 



$„(M^mi,m 2 ,...,m„) = J dM^_i#2(M^, M n _i, m„)$ n -i(M^_i, mi, mi, . . . , m n -i) 



(19) 



,mi,m 2 , . . . ,m„-i) 



4M n 



(M„-m„) 



dMl_ x 



\(Ml,Ml_ x ,ml) 



8M2 



, mi, m2, ■ • • , rn n -i ) 



(E^i 1 ™;) 2 



with the integration limits on M„_i following from its definition in Eg. 1151 It has to be emphasized that the angles in dfi* are 
each time calculated in the centre-of-mass system of ki with the invariant mass Mi. The resulting recursion relation is clearly 
of advantage when describing cascade decays of particles k n — > k n ~ip n — ► fc n -2,Pn,Pn-i — > . . .; it also proves that the 
n-particle phase space of Eq.|5|can be reduced into a sequence of two-particle phase space terms, as shown in Figure|3| 



Pn Pn-l Pn-2 



Pi+1 Pi 




Fig. 3: The diagrammatic representation of consecutive s-splits. 
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It can further prove of advantage to loosen up the splitting terms of Eqns. H6ll7l so that instead of summing to n-1 one 
groups an arbitrary set of £ particles: 

1 = J dM?5(k? - M?)&(k?), (20) 

/n 
d 4 fc ; 5 4 (p-fc ; - pj). ( 21 ) 
i=i+i 

which, when repeating the procedure in recursion relation of Eq. fT9l results in an expression: 

(A/ ! + 1 -m i + 1 ) 2 

$„(M^,mi,m 2 ,...,m„) = / dM, 2 $„_i + i(M^, Mi, mj+i, . . . , m n )$i (Mi, mi, m 2 , . . . , mi), (22) 

and thus effectively splits the phase space into two subsets, equivalent to introducing an intermediate(virtual) particle with 
momentum fc;. 

The number of splitting relations and the number of particles in each group as given in Eq. 12 II can be chosen in any 
possible sequence, thus meaning that the grouping sequence is arbitrary and can be adjusted to fit the topology in question. 3 

At this point some modifications were introduced to the algorithm in order to adapt it to the specifics of the processes 
expected at the LHC. Kajantie and Byckling namely assumed that the generation sequence would be 'down' the cascade (i.e. 
by sampling first a M n value, then M n _i value etc. as is indeed most often done in Monte-Carlo Generators). This might 
however not be optimal in the LHC environment since the available centre-of-mass energy for the hard process (s) can vary 
in a wide range of values (c.f. Equation |2j and has to be sampled from a distribution itself. The shape of the distribution 
function for s is expected to behave as a convolution of the peaking behavior of all participating invariant masses times the 
parton density functions (c.f. Eq. it subsequently seems to be more natural (and efficient) first to sample the individual 
propagator peaks and then their subsequent convolutions. Furthermore, by generating the invariant masses 'up' the cascade 
(i.e. first M2, M3 .M n and finally s) the kinematic limits on the branchings occur in a more efficient way (bound on the \/A 
values, see Eauations l 14l and l24l . which is very convenient since in the LHC environment no stringent generation cuts should 
be made on the inherently non-measurable s as it cannot be accounted for by an analogous cut in a physics analysis. 

A necessary modification of the algorithm would thus be to reverse the generation steps by starting with the last pair(s) of 
particles. In terms of integration (i.e. sampling) limits this translates into changing the limits of Eg. 1191 

$„(M 2 ,mi,m 2 , . . . ,m„) = (23) 

U m " )2 JMMlM^ml) , 

dM -- — m l dUn 

(EIL - ! 1 ™.) 2 

(M„_i-m n _i) 2 



JxiM^M^ml,,) 

dMl ->- smT, / '«-'» 1 

(rr=i 2 ™.) 2 



dMU± — / ,/<»; 

(E}=i™ 3 ) 2 

(M 3 -m 3 ) 2 



dM WHMlMJ, m $ 



12 VI/-* 

(mi+m 2 ) 2 



8M| 



■~J\{Ml,m\,ml) 
8M| 



3 Suggestions of 1151 on how to pick random number sequences will not be used since one might like to couple this method with an adaptive algorithm to 
improve the sampling efficiencies. 
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which accommodates the mass generation sequence: k n — > fc„_i +p n — ► .. . (i.e. first sample M„_ 1 , then M n _2 etc.), into 



(M„-m n y 



Jx(MlM^,ml) 

$ n (Mlm u m 2 ,...,m n ) = J dM 2 n ^ — / d£l* n (24) 

(M„_ 2 +m 7l _ 1 ) 2 
(M n —m n — m n _ 1 ) 2 



2 ^(Mj^Mj^m 



2 ) 



(M„_ 3 + m„_ 2 ) 2 

(Afn-Ej= i+ i m 3 ) : 



(JVIi-l+m;) 2 
(M»-E" =3 ™i) 2 



/ 



8M| 

where one first samples the mass M2, M3.M n _i in the appropriate limits. 

In some topologies symmetric cases of mass generation can appear (as shown in Figure^} where the integration sequence 
is ambivalent (e.g. in Figure Q] the ambivalence is which top quark invariant mass to generate first.) and after a choice is 
made (since one of the two cases in the symmetric pair has to take precedence) the procedure itself remains not entirely 
symmetric. Detailed studies have shown that it proves useful to include all permutations of such ambiguous sequences into 
the MC algorithm in order to 'symmetrise' the solution and thus make it easier to process by further additions (e.g. adaptive 
algorithms). 

2.2 The t-type Branching Algorithms 




Fig. 4: A diagram of a generic 2^2 t-channel process. 



The t-splits are a specialty of the KB formalism due to the advanced calculation of the limits on the (massive) t-variable. 
The formalism can be introduced by observing that in case of a p a + pb — > pi + P2 scattering the momentum transfer is 
characterised by the (Mandelstam) variable t = (pi — p a ) 2 (c.f. Fig|3J. It is thus sensible to replace the dfll — dcos 6* dtp* 
integration in the two body phase space integral of Eq. El with integration over the t variable. Writing the definition of t in 
the centre-of-mass system one gets: 

t = q 2 = (Pa- P i) 2 (25) 
= ml + rn{ - 2E* a E{ + 2p* a p\ cos 6{ 

and hence: 

dt — 2p* a p\d cos 6* (26) 
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Using the latter substitution together with Eq. I12ll3l and the analogue for p 



Pa = 



one obtains in place of Eq. 1121 



$2(s,mi, ma) 



y , A(s,mg,mg) 
2^ 

Pi(s,mi,m 2 ) 



(27) 



(28) 



dt dip* 



A-^f \{s,ml,ml) 



+ 2lT 

dt / dip* 



With the integration variable change the integration domain changes from [—1, 1] fordcos#* to [t ,t + ] for the dt integration. 
The limits are obtained by inserting the cos 6* limits into Equation l26l 



=ml + ml- 2E* a E{ ± 2p* a p* 



or in the Lorentz invariant form (c.f. Eq. 112113 

,± 



(29) 



(30) 



2 . 2 (s + m a - m b )(s + mi - ra 2 ) 

m a +m 1 - ^ ^ i 

2s 

V / A(s,mg,mg)A(s,mj,m|) 
2s 

As a step towards generalisation one has to note that the kinematic limits t can also be derived from the basic four- 



particle kinematic function G(x. y.z.u. v. w) \ 1611151 . where the function G can be expressed as a Cayley determinant: 



G(x,y,z,u,v,w) = -- 
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The kinematic limits on t are in this case given by the condition 



G(s,t,ml,m 2 a ,ml,ml) < 0, 



(31) 



(32) 



it should be noted that the above condition gives either t limits given a fixed value of s or equivalently s limits given a 
fixed t value. 

In search of a recursion relation involving t-variables one can note that in Eq. I19l the angle in cos 8„ is equivalent to the 
scattering angle in the centre-of-mass system of the reaction p a + Pb — > k n _i + p n and thus given by: 



t n 



with the t n _ 1 limits expressed by: 



= (Pa — fcn-l) 

= ml + M^_i - 2Elk°-i + 2p^fc*_i cos0£_i 



G(M n ,t n - 1 ,m n ,m a ,m b ,M n _ 1 ) < 0, 



(33) 



(34) 



and the p*. given by Eq. 1271 In order to produce a more general picture it can further be deduced that the next angle in the 
recursion 8* l _ 1 , is the scattering angle of the subsequent process p a + (pb — p n ) — > k n _2 + p n -i in the centre-of-mass 
system of fc„-i; the (pb — p n ) = qn-i is in this case considered as a virtual incoming particle with momentum q n -i (c.f. 
Figure |5j. 

It immediately follows that for a general process p a + — > ki + Pi+i with: 



qi — Pb- 22 P] = Pa ~ ki ' * = tl ' In =t n = m b 



(35) 
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q;+i 



k, 



Fig. 5: The diagrammatic representation of the method applied in translating the multi-(peri)pheral splits 
into a 2 — > 2 t-channel configuration. 



a general expression for tj becomes in the centre-of-mass frame of ki+i: 

ii = {.Pa fci) 



m 2 a + Mf - 2E* a 



,(i + l),o»{i+i) + 2 p* a ^^k'^ ± > cose* 



where momenta in centre-of-mass frame of k;+i, denoted with the superscript *(i + 1), are given by: 



A: 



:(i+t) 



*(<+!) 



^\{Mf +1 ,Mf,m} +1 ) 
2M !+ i 



(36) 



(37) 



(38) 



and the corresponding energies fc^ 1 * and Ea can simply be obtained by using the analogues of Eauations ll2ll3l or the 
usual Einstein mass-energy relations directly. The corresponding tj limits given by: 

G(Mi +1 ,U,m 2 i+1 ,ml,t l+1 , M?) < 0, 

and the recursion relation of Eq. ll9l becomes: 

$„(M^mi,m 2 ,...,m„) = 

(M„-m„) 2 ^ 2vr 'n-l 

: / d^i* / dt„-i $ n _i(M^_i,mi,m2, ...,m„_i) 



(39) 



(40) 



(Efei 1 ™.) 2 



As already argued the resulting set of (sj = Mf,ti) can again be sampled in any direction with respect to the cascade 
by applying the appropriate change in the integration limits (c.f. Eq. 1191 and l24l . The recommended approach (i.e. the 
introduced modification of the algorithm) is again to first sample the invariant masses in the reverse cascade direction (i.e. in 
the sequence M2, M3, . . . , M n ) and then the i, values within the limits calculated from Eq. B4l down the cascade (i.e. in the 
order of t n _i,t n _ 2 , . . . ,ti). 

To sum up, it has been shown that using the Kajantie-Byckling formalism the phase space for any topology can be split 
in a set of s-type and t-type 2^2 branching steps (modules) given by recursive formulae of Eauations l24l and l40l 
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3 Propagator Sampling 



A well known theoretical issue is that one can expect the most prominent peaks in the differential cross-section of a specific 
process in the phase space regions of high propagator values in the corresponding probability density. Consequently, in the 
scope of complementing the modular structure of the derived Kajantie-Byckling based phase space sampling, new approaches 
were also developed concerning the numerical sampling methods of the relevant kinematic quantities. 

In order to get small variance in the Monte Carlo procedure one would thus like to include the appropriate peaking 
dependence of the relevant momentum transfers q 2 in the importance sampling function. It however turns out that since the 
momenta transfers q participate also in the propagator numerators (typically in p M q M /q 2 ) and since in process of interest 
one mostly finds several Feynman diagrams contributing to the final probability density, thus causing interferences, it is very 
difficult or even impossible to estimate the exact power of momenta transfers in the sampling functions for different propagator 
peaks. In other words, the probability density dependence on the momentum transfer q 2 can in general be approximated with 
the dependence 1 / (q 2 )" where the best value of v must be determined separately (on a process by process basis). 

In view of the latter, general formulae have been developed for sampling the x~ v shape II 2111 71 : Given a pseudo random 
number r G [0, 1] and limits x G [x~, x+] the value x distributed as x~" is obtained from the formulae in Eq. [8]as: 

x = [xZ u+1 -(l-r)+xZ u+l -r]~ 1 fc; v =/= i ; (41) 
x = ^± T ; v = \. (42) 

Using the analogous (unitary) approach a recipe for resonant (Breit-Wigner) propagator contributions of the type: 

BW ® = (s-m4 + M^ (43) 
with s G [s_ , s+] and a pseudo random number r G [0, 1] is available by the prescription: 

s = M 2 + MY • tan [(u+ -U-)-r + u-] (44) 

« = atan I Mr J (45) 

Following similar arguments as for the non-resonant propagators one can surmise that the best sampling function for 
resonant propagators could in general be a Breit-Wigner shape modified by a factor s", v G [0, 1]. In 1 10 1 it was found that a 
shape: 

BW ^ = FWT1¥ (46) 

works quite well for a set of processes and a corresponding sampling recipe was developed. In addition, studies in 1 18 show 
that a resonant -^sx Breit-Wigner shape: 

BW ^ = (s-M^ + MV* (4?) 

should be expected in a range of decay processes. Detailed studies have shown that it is in general better to introduce a 
s" ,v G [0, 1] dependence even if it over-compensates the high mass tails of the corresponding differential cross-section 
distribution since this provides an overall reduction of the maximal weight fluctuations in the Monte-Carlo event generation 
procedure. 



3.1 The Inclusion of Mass Effects in Propagator Sampling 

Studies have shown that the x~ v approximation works quite well for t-channel type propagators since the phase space suppres- 
sion factor \f \ participates in the denominator, as shown in Eg. 1401 and thus contributes only to the x~ v slope. Contrariwise, 
while the x~ v approximation still works reasonably well when describing the s-channel type propagators involving particles 
with high virtuality and/or decay products with low masses, it can be shown that this is not necessarily the case in the LHC 
environment, where the presence of massive decay products can significantly affect the invariant mass distributions. As it can 
be seen in Figure|S|the shape of the propagator dependence can be strongly suppressed by the phase space \f\ (boost) factor 
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at low values; thus the sampling function approximation for non-resonant propagators could be approximated with something 
like: 



r r \ _ y --\"! -ai ■■•bl J- _ V -ai ,.qs 

./nr(s) ; —; ^zr ( 48 ) 



^/A(s,m 2 ,m 2 ) 1 _ ^/\(s,mj. mr 



and similarly 



= V / A(a,mg,mg) y/s = y , A(s,mg,mg) 

M[ ' s ' (s-M 2 ) 2 + M 2 r 2 ^i-((s-M 2 ) 2 + Af 2 r 2 ) 1 ; 



for resonant propagators. 

As it turns out the two functions cannot be sampled by the well known unitary algorithms (i.e. the biggest collection of 
recipes 1 17 1 yielded no results); already the integral values of the functions yield complicated expressions which cannot be 
easily calculated, let alone inverted analytically. The solution was to code numerical algorithms to calculate the integrals (i.e. 
cumulants) explicitly. 

After the integrals are calculated, their inverse and the subsequent sampling value can again be obtained numerically. 
Namely, resorting to the original definition of the unitarity sampling recipe in Eq.[8] by replacing the normalised gi(x) with: 

- J {X) ■ (50) 
/ f(x)dx 

which in turn gives: 

x x + 

J f(x)dx = r- J f(x)dx, (51) 

X _ x _ 

where f(x) is the non-negative function one wants to sample from, [x_,x+] is the range of values of the parameter x we 
want to sample and r a pseudo random number r G [0, 1]. As already stated (c.f. Eq|8j, in the case when the integral of 
the function f(x) is an analytic function, F{x) = f(x) dx, and has a known inverse F _1 (x) one can construct explicit 
unitary prescriptions by: 

x = F- 1 (r-[F(x+)-F{x-)] + F(x-)) (52) 

as given for two particular cases in Eq. 1421451 

In the cases the integral can not be inverted, the prescription of the Eq. 151 l ean directly be transformed into a zero-finding 
request; thus, since both the integral and the first derivative (i.e. the sampling function and its cumulant) are known, the 
Newton-Rhapson method is chosen as the optimal one for root finding: 



1x - + 

J f(x)dx - r ■ J f(x)dx i = 
X _ X _ 

(x x + 

J f(x)dx -r- J ] 



(53) 



With a sensible choice of starting points the procedure generally takes on the order of ten cycles until finding the root with 
adequate numerical precision. The overall generation speed is still deemed quite reasonable. 
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The integration of the phase-space suppressed resonant propagator of Eq. l49l vields a rather non-trivial expression: 

/ R (s) ds = 



\/X(s,ml,ml)ds 



(m a +m b ) 2 



(m a +m(,) 2 



/i-((s-M 2 ) 2 + M 2 r 2 ) 



(55) 



y (s — a)(s — b) ds 
yfs- ((s-M 2 ) 2 + M 2 T 2 ) 



-2ia6T 



-&TM 2 (T 2 + M 2 ) 

iarcsmh( — — ), — 
Wa b 



i arcsinh(- 



—6. a 

7T ], b 



+ (i T + M) (a + i (r + i M) M) (b + i(T + i M) M) II 
+ (V + i M) (b + (— i T — M) M) (ia + (T — i M) M) II 

- (i r + M) (a + i (r + i M) M) (b + i (r + % M) M) n 

- (T + i M) (b + {-iY - M)M){ia + {T -i M) M) U 

where the variables a, b stand for a = (m a + nib) 2 and b = (m a — nib) 2 and the functions F[tp, k] and H[tp, k, n] are 
the Legendre's incomplete elliptic integrals of the second and third kind with complex arguments. In order to perform the 
calculations the latter functions had to be coded from scratch since they were not found in any (publicly available) computer 
libraries or code repositories. The prescriptions for calculating them were found in 1 19 1; the results were checked against the 
values given by Mathematica™. 

In the special case m a = nib the above expression simplifies into: 



M (-i T + M) 
b 


i arcsinh(- 


J-b 
\/a 


' 6 


~M(iT + M) . 


arcsinh( — 


^) 


a 


[ b 


fa h 


b 


M{-iT + M) 
b 


i arcsinh(- 


J-b 
V~s 


' b 


' M (iF + M) . 
b 


arcsinh( — 


^) 


a 
b 



f , , s , f \/Ms,m 2 a ,m 2 a )ds 

J Ms)dS = J ^.(( S -M 2 ) 2 +M 2 



(2m a ) 2 



(2m a ) 2 



r 2 ) 



(56) 



\J (s — a) ds 



•((s — M 2 ) 2 + M 2 F 2 ) 

a 

l 

TM s/a + i-iT - M) M 
(i a + (T - i M) M) arctan( 



Ja + (-iT - M) M 



i a/o + (-iT- M) M Ja + i (T + iM) M arctan(- 



V— a + z 



y/a + i (T + iM) M 

The result of integrating the phase-space suppressed non-resonant propagator (Eg. 1481 yields a similarly non-trivial result: 



(ma,+m b ) 



J /nr(s) ds = J 

(m a +m b ) 2 



\A(s,m 2 ,m 2 ) ds 



(57) 



-2 y/(a -s)(b- s) Fx [-„, - (|) , - (I) , 1 - u, | , f ] 



+ 



^^/(-a + b) (a-a)T[l-u] F [-»/,-(§), f -«/, f] 

o-yr^rfl-i/] 
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where the function F[a, f3, 7, x] is the Gauss Hypergeometric function and the Fi [a, j3, /3' , 7, x, y] is the two-parameter 
(Appell) Hypergeometric function |20|. Both functions can be calculated by using the prescriptions in |20|; it however 
turns out that the calculation of the Fi [a, (3, /?', 7, x, y] to a certain (high) precision is almost two times slower than the 
explicit numerical calculation of the integral to the same precision. Subsequently the numerical evaluation of the Gauss 
Hypergeometric function F[a, /3, 7, x] was retained since it participates in the m a = nib simplification and the calculation of 
the integral was done by using a 50-point Gauss-Legendre quadrature with weight function; the weights were calculated 
by 1211 . The implementation of the (Appell) Hypergeometric function calculation was used as a cross-check to confirm the 
correct implementation and precision of the numerical method. 
As already mentioned, the above integral again simplifies for m a = nit,: 



s 



s 




(58) 



S 




a 



2 





1 5 



and the Gauss Hypergeometric function F[a, j3, 7, x] is in this case calculated by the methods described in |20| with some 
improvements analogous to the ones described e.g. in 1221 . 
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4 Application of the Method 

The phase space 'modularisation' described in this paper has successfully been applied for phase-space generation in the 
AcerMC 2.0 Monte-Carlo generator! H ]■ The program uses the multi-channel phase space generation where each channel 
corresponds to an expected phase space topology as derived from the participating Feynman diagrams. In the AcerMC 2.0 
this information was obtained from the modified MadGraph|6| program which also supplied the probability amplitudes for 
the implemented processes. Each channel topology was in turn constructed from the t-type and s-type modules and sampling 
functions described in this paper together with some additional importance sampling techniques for space angles and rapidity 
distributions described in detail elsewhere |23 12l ll3lfTol . The unknown slope parameters (denoted v in the text, c.f. Eq. l48> 
of the invariant mass sampling functions for non-resonant propagators were obtained by short training runs of the program on 
a process by process basis. 

As a further step the multi-channel self-optimisation procedure was implemented in order to minimise the variance of the 
event weights further |7|. A few representative invariant mass distribution comparisons between the implemented sampling 
functions and the actual differential distributions are shown in Figure|6| 



■ Sampling 

■ da/dm. 




^ Sampling 




- — da/dm wbb 






> 



in 20 30 40 50 60 70 80 90 100 
m tb [GeV/c 2 ] 




m,, [ GeV/ c- ] 



Fig. 6: A few representative invariant mass distribution comparisons between the (normalised) sampling 
functions and the normalised differential cross-section as obtained with AcerMC 2.0 Monte- 
Carlo generator. Left: The invariant mass of the bb pair in the process ud — ► W + g* — > l + ^ibb. 
Center: The invariant mass of the Wbb system (equivalently the hard centre-of-mass energy %/§) 
for the same process. Right: The invariant mass of the 11 pair in the process gg — > Z° /"f*bb — > 
Ubb. All the distributions were obtained using the prescriptions of this paper without the adaptive 
algorithms also used in the AcerMC 2.0 Monte-Carlo generator. As one can see the approxima- 
tions used seem to work quite well. 



As a further estimate of the success of the methods a comparison of the variance in the differential cross-section deter- 
mination are presented in Table[Jfor AcerMC 2.0 | ll] and an earlier version, AcerMC 1 .4 1101 . which uses the standard 
phase space sampling techniques 1 121 1 1 31 l24l 1251 1 101 . The comparison is done for a few representative processes. It has to 
be stressed that in AcerMC 1 .X versions the phase space was constructed by hand on a process by process basis meticulously 
tuned while in the new AcerMC 2.0 the 'automated' approach sufficed. A further improvement was the inclusion of the 
ac-VEGAS algorithm 1 10 1 which reduces the maximal event weights in order to improve the unweighing efficiencies; the 
comparison between the unweighing efficiencies reached in the standard and new approaches (AcerMC versions 1 .4 and 2.0) 
is also given in Table[J The event weight variance V a presented in Tablets for N measurements of the weights wi defined 
as: 



V " = N±l (59) 



Likewise, the customary definition is used for the unweighing efficiency estimate: 

N 



(60) 

^max 
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where the w max is the maximal event weight obtained in N trials. It also might be relevant to stress that since the average 
weight is the total cross-section estimate for the considered process a ~ 53i=i W i/N the only quantity that is allowed to 
change in order to improve the unweighing efficiency is the maximal weight w max . 

Table 1: The process cross-section variances with their uncertainties and unweighing efficiencies obtained 
for a few implemented processes basing on the old and new phase space sampling techniques 
as implemented in AcerMC 1 .4 and AcerMC 2.0 respectively. The variances are given for 
a sequence of 10 weighted events (i.e. algorithm iterations) obtained by using the procedure 
described in the text. The unweighing efficiencies were estimated from samples containing ~ 10 6 
weighted events. 



Process 




AcerMC 2.0 V CT [pb 2 ] 


AcerMC 1 .4 V CT [pb 2 ] 


AcerMC 2.0 e 


AcerMC 1 .4 e 


99 - Z/(-» 


U)bb 


0.129 ■ 10~ 2 ± 0.52 ■ 10" 5 


0.159- 10" 2 ±0.61- 10" 5 


37% 


33% 


qq^ W(-> 


iv)bb 


0.390 -10~ 2 ± 0.15 • 10 -4 


0.533- 10~ 2 ±0.18- 10 -4 


35% 


33% 


gg —* tibb 




0.522 ■ 10" 4 ± 0.19 • 10" 6 


0.972- 10~ 4 ±0.44- 10" 6 


36% 


20% 



Table 2: The process cross-sections and variances with their uncertainties and unweighing efficiencies as 
obtained for two sample 2^6 processes implemented in AcerMC 2.0 Monte-Carlo generator. 
The results show that the sampling procedure scales quite efficiently with the increase in the 
number of Feynman diagrams and sampling channels. Detailed studies show that with inclusion 
of the new propagator sampling the angular distributions become the restraining factor. The 
cross-sections and variances are given for a sequence of 2 ■ 10 5 weighted events (i.e. algorithm 
iterations) obtained by using the procedure described in the text. The unweighing efficiencies 
were estimated from samples containing ~ 10 6 weighted events. 



AcerMC 2.0 Process 




a[pb] 


V CT [pb 2 ] 


e 


gg^ti^ bbW + W~ -> Mutlut 


(3 Feyn. diag./2 sampl. chan.) 


4.49 


0.80- 10" 4 ±0.39- 10" 6 


14% 


gg -> bbW + W- -s- bblvtivt 


(31 Feyn. diag./13 sampl. chan.) 


4.77 


0.77- 10~ 4 ±0.29- 10" 5 


17% 



The improved and automated phase space handling provided the means to include the 2^6 processes like e.g. gg —* 
ti — > bbW + W~ — ^ bblviqxqy. (c.f. Fig. [3 which would with the very complicated phase space topologies prove to be too 
much work to be handled manually. The studies show that the overall unweighing efficiency which can be reached in the 
2^6 processes by using the recommended phase space structuring is on the order of 10 percent. As an example, process 
cross-sections and variances with their uncertainties and unweighing efficiencies as obtained for two sample 2^6 processes 
implemented in AcerMC 2.0 Monte-Carlo generator are presented in Tableland the corresponding weight distributions are 
shown in Figure|7| The two processes share the same initial and final state but include different subsets of Feynman diagrams; 
the gg — > ti — > bbW + W~ — > bblviivi process includes only the three diagrams containing the ti intermediate state while in 
the second process gg — > bbW + W~ — > bblutlut all 31 diagrams involving the bbW + W~ intermediate state are included. 
Subsequently, two sampling channels were constructed to describe the topologies in the first and 13 sampling channels in the 
second process of the two described above. The results show that the sampling procedure scales quite efficiently with the 
increase in the number of Feynman diagrams and topologies. Detailed studies show that with inclusion of the new propagator 
modeling described in this paper the sampling of angular distributions becomes a new restraining factor in achieving optimal 
weight variances and unweighing efficiencies which might therefore be worth investigating further. 
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Wt x 2 • 10 6 [ pb] 



Fig. 7: The weight distribution of the sampled events for the gg — > tt — > bbW + W~ — > bUviivi (light 
gray histogram) and gg — » bbW + W~ — * bblvilvt (black histogram) processes as obtained with 
AcerMC 2.0 Monte-Carlo generator. One can observe the well defined weight range for the two 
processes; as it turns out the weight distribution is even marginally better for the (more complex) 
second process, possibly because the higher number of sampling channels manage to cover the 
event topologies in phase space to a better extent. 



5 Conclusion 

In this paper the revised Kajantie-Byckling approach and some improved phase space sampling techniques for the massive 
multi-particle final states were presented. In order to adapt the procedure to the LHC environment the modifications necessary 
for reversing the sampling order were introduced and new invariant mass sampling methods, which attempt to describe the 
propagator dependence of the probability density together with the phase space suppression due to the presence of massive 
particles, were developed. 

The developed procedures have been implemented in the AcerMC 2.0 Monte-Carlo generator! 1 1 1- Basing on the encour- 
aging evidence provided by the AcerMC 2.0 implementation of the approach it seems reasonable to argue that the methods 
presented in this paper should substantially simplify and automatise the phase space integration (sampling) techniques while 
retaining a respectable weight variance reduction and unweighing efficiencies provided by the most advanced phase space 
sampling techniques developed so far I12I H3 24 25 10 J. Furthermore, the authors believe that these techniques could easily 
be combined with algorithms of the type Sarge 1 4 1 or HAAG 1 5 1 to provide successful sampling of the (~massless) final state 
particles as e.g. final state gluon radiation. 
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